IASSNS-HEP-94/80; cond-mat:9410007 
A Cellular Automaton Model for Diffusive and Dissipative Systems 



T. C. Chan\ H. F. Chau^'^*, and K. S. Cheng^ 
^ Department of Physics, University of Hong Kong, Pokfulam Road, Hong Kong. 
^ School of Natural Sciences, Institute for Advanced Study, Olden Lane, Princeton, NJ 08540, USA. 
^ Department of Physics, University of Illinois, 1110 West Green Street, Urbana, IL 61801, USA. 

(February 1, 2008) 



Abstract 



We study a cellular automaton model, which allows diffusion of energy (or equiv- 
alently any other physical quantities such as mass of a particular compound) at 
every lattice site after each timestep. Unit amount of energy is randomly added 
onto a site. Whenever the local energy content of a site reaches a fixed threshold 
Eel, energy will be dissipated. Dissipation of energy propagates to the neighbor- 
ing sites provided that the energy contents of those sites are greater than or equal 
to another fixed threshold £'c2(< Ed)- Under such dynamics, the system evolves 
into three different types of states depending on the values of Ed and Ec2 as re- 
flected in their dissipation size distributions, namely: localized peaks, power laws, 
or exponential laws. This model is able to describe the behaviors of various phys- 
ical systems including the statistics of burst sizes and burst rates in type-I X-ray 
bursters. Comparisons between our model and the famous forest-fire model (FFM) 
are made. 
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I. INTRODUCTION 



For many years, cellular automaton models have been playing important roles in the study of 
some non-linear and spatially extended systems. In most cases, these models cannot be solved ana- 
lytically and their investigations rely mainly on computer simulations. Forest-fire model (FFM) [Q, 
Eden model and earthquake simulations are some examples. Besides, a lot of simplifications 
has been taken in these models, making their underlying physics unclear. In this paper, we intro- 
duce a simple cellular automaton model that accounts for the stochastic energy (or equivalently 
any other physical quantity) introduction, occasional energy dissipation, and energy diffusion in an 
open physical system. This kind of non-equilibrium systems are common in nature. We find that 
the system shows different behaviors with different parameters. And with an appropriate choice of 
parameters, the model can be applied (sometimes after a small changes in the geometry) to explain 
the statistics of the occasional outburst of energy (or other physical quantities) in various physical 
systems including the thermonuclear runaways in X-ray bursters |^ and the sudden CO2 gas release 
in crater lakes |^. 

We introduce the model in Section II. Then we discuss the significance of various parameters 
and present the general results of computer simulation in Section III. In Section IV, we concentrate 
on a special case and compare it with the FFM. Finally, a conclusion together with a discussion on 
the application of our model can be found in Section IV. 



We consider a L x L square lattice each of size AL x AL. The energy contained in a lattice site 
r is denoted by E{r), which is a non-negative real number. Since the sites are of equal areas, the 
ratio of energy to energy density is a constant for all site. Evolution of the system is governed by 
the following rules: 

1. Energy Introduction: at each timestep, a unit amount of energy is added to a randomly 
selected site tq, i.e. 



It should be noted that the energy introduction rate to the entire system is being fixed. Our 
energy introduction rule takes only spatial fluctuation into account. 

2. Triggering Of Energy Dissipation: whenever the energy of a site fq reaches a fixed value E^i 
called the triggering threshold, energy stored at that site will be dissipated, i.e. 



Note that the triggering threshold Ed is the same at all sites and is time independent. The dis- 
sipation process is called "burning" . Since all sites are of equal areas, we can regard the "fire" 
as being triggered whenever the local energy density at some place exceeds a predetermined 
triggering threshold. 



II. THE MODEL OF DIFFUSION AND DISSIPATION 




(1) 




(2) 
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3. Propagation Of Energy Dissipation: nearest neighbors of a burning site burn provided that 
their energy content are greater than or equal to a fixed value Ec2 (< E^i) called the propa- 
gation threshold, i.e. 



if E{r'o) > Ec2 ,then E{t'o) 



(3) 



whenever Fq is a nearest neighbor of the burning site fq. We use periodic boundary conditions 
in the determination of neighboring sites. One may regard the triggering threshold E^i as the 
(energy) density required to activate the dissipative process. Once the "burning" takes place, 
the energy release may heat up the neighboring sites and hence the density threshold required 
for the fire to propagate, Ec2, is reduced. This is valid in many physical and chemical reactions. 

The burning process continues until no more lattice site in the system catches fire. The energy 
dissipation (or outburst) is then completed. This is the fastest process in the system and so 
the whole burning process is assumed to take place in a single timestep. 

Finally, we define a "burnable cluster" as a collection of sites such that all of them will catch 
fire in case any single one of them burns. This is a useful concept when discussing the statistics 
of energy outburst. 

4. Energy Diffusion: diffusion takes place in each timestep, i.e. 



where AE depends on the system configuration, the diffusion constant D, the lattice size AL, 
and the physical time At corresponding to a cellular automaton timestep. But without lost 
of generality, we can always rescale our time and space to make both At and AL equal 1. 
Furthermore, the method we use to evaluate AE{r) is reported in the Appendix. 

Various initial conditions, such as starting with all E{r) = 0, have been used. However, the 
long term statistics is independent of the initial conditions because both burning and diffusion erase 
the history of the system. Although the model we have introduced is based on a two-dimensional 
square lattice, generalizations to other dimensions and griding methods are straight forward. More 
complete study of the model in various dimensions will be reported in future works 0. Finally, a 
snapshot of a typical system configuration immediately before a fire is triggered is shown in Fig. |I[ 



Before presenting the results of our numerical simulation, let us try to argue qualitatively what 
we expect to see. 



There are three tuning parameters in our model: namely Ed, Ec2, and D. Large average 
dissipation size is a direct consequence of having large burnable clusters immediately before the 
system catches fire. And large burnable clusters can be obtained in either one of the following 
ways: 




for all r 



(4) 



III. EXPECTATIONS AND SIMULATION RESULTS 



A. Expectations 
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1. Large Diffusion Constant: diffusion spreads energy out and wipes out information on where 
the energy packet is first introduced to the system. Having a large D means that immediately 
before a sudden energy dissipation, the energy content in each site is approximately equal 
except for those few sites where energy packets are recently introduced. So, we expect the 
whole system to be covered by a single burnable cluster immediately before each burning. 

2. Small Propagation To Triggering Thresholds Ratio: in this case, most of the sites have enough 
energy content to continue the propagation once the burning begins. Once again, we expect 
to find a single burnable cluster covering the entire system immediately before each burning 

i- 

Using the same argument, small average dissipation size can be obtained when diffusion is 
unimportant and when Ec2 is comparable to E^i- For every fixed and Ec2, the average burnable 
cluster size immediately before a burning increases with D. Therefore, starting from a subcritical 
system (i.e. system with small burnable clusters only) with D = 0, a supercritical system (i.e. 
system with large burnable clusters only) can be obtained simply by increasing the value of D. We 
also expect to find a critical value of D, which is a function of E^ and Ec2, such that burnable 
clusters of all sizes can be found immediately before burning. This is the critical state in our model. 
Numerical experiments we have done confirm our expectations, and we are going to report it in the 
coming Subsection. 



B. The roles of diffusion 

In order to study the effects of diffusion, we vary the diffusion constant D from 10~^ to 10~^ in 
a 64 X 64 lattice with E^i and Ec2 equal 5.0 and 2.0 respectively. We found that the average size 
of energy dissipation (S) increases as the diffusion constant D increases (see Fig. |]). Although we 
only report the system behavior when it assumes the above parameters, the general behavior of the 
system using different values of thresholds are unchanged. 

We observe that when D ^ 3 x 10~^, both the energy dissipation size S and the time since last 
dissipation T become very regular (see Fig. ^) . This is the consequence of having a single burnable 
cluster covering the whole surface just before the burning. Thus all the energy available in the 
system is dissipated during the burning. As a result, S equals the total amount of energy dumped 
into the system since its last dissipation which is in turn equal to T. The strong correlation between 
S and T for high diffusion system is shown in Fig. ^ (the circle symbols) . To study this correlation, 
we define the correlation coefficient F by: 

r^m^mii, (5) 

aso-T 

where as and ctt are the standard deviations of S and T respectively. It can be shown that F 
varies from -1 to 1, and the value of 1 (-1) implies that S and T are completely positive (negative) 
correlated while means that they are uncorrelated. The dashed line in Fig. |^ shows that F increases 
with increasing D indicating that a strong positive correlation between 5* and T for system with 
large diffusion constant. In this case, the distributions of energy dissipation P{S) and time interval 
P(T) are the same, and both of them show localized skewed peaks (see Fig. 
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For a smaller D {5 x 10~^ ^ D ^ 3 x 10~^), the strong correlation between 5* and T begins to 
break down. Weaker diffusion permits greater ffuctuation in the energy content at different sites. 
The system is likely to catch fire before all the sites are connected. Hence not all energy available 
in the system is dissipated, and some clusters may leave after a burning. As shown in Fig. ^, {S, T) 
[the plus symbols] are distributed around their mean values. In addition, if no burning is triggered 
for sufficient long time, a system-wide dissipation can happen. F decreases from 1 in this region as 
well. A consequence of the incomplete burning is that the probability of having a triggering site 
in a small burnable cluster is greatly increased. Events with small energy dissipation size begin to 
appear and hence the correlation F increases slightly as D decreases around D = Q x 10~^. 

As D decreases further (10~^ ^ D ^ 5 x 10^^), small dissipation begins to dominate. In fact, 
P{S) changes from localized peak to power law and finally to exponential decay. In our simulation, 
we found that the critical diffusion constant Dcru to be 1.4 x 10~^. Surely this critical value is 
a function of Ed and Ec2- For even smaller D {D < 10~^), diffusion becomes insignificant. The 
heights of nearby sites are nearly uncorrelated and hence P{S) shows a very early cutoff (i.e. an 
exponential decay). Fig. ^ depicts the localized skewed peak, power law and exponential behaviors 
of P{S) as a result of different diffusion constant D. 

It should be noted that when P{S) follows a power law, the corresponding distribution of time 
interval between successive burning P{T) decreases exponentially for large T. It means that there is 
a characteristic time interval for the occurrence of burnings even when scaling behavior is observed 
in the distribution of dissipation sizes. It is because the burnable clusters are spatially separated 
and the largest burnable cluster is only a portion of the system. The total number of critical sites 
in the entire system and hence the probability of triggering a burning are more or less the same all 
the time. It is expected that this case has similar features with respect to the one studied in the 
next section. 

C. The roles of triggering and propagation thresholds 

In our simulations, we have not studied the effect of the thresholds on the model in detail. 
However, we believe that the above behaviors would be also observed when we vary Ed or Ec2 
instead of D. In fact, when one of the three parameters is fixed, we may observe any one of the 
above behaviors by carefully choosing the other two parameters together with a reasonable lattice 
size. For example, reducing Ec2 has similar effects as increasing D (see Fig. |^). 

IV. COMPARISON WITH FOREST-FIRE MODEL 

Forest-fire model (FFM) was first introduced by Bak et al. as an example of self-organized 
criticality (SOC) [0]. They considered a lattice of cells where individual cell belongs to either one 
of the following states: a green tree, a burning tree, and a void (i.e. no tree). In each timestep, a 
tree grows on an empty cell with probability p, and a green tree burns if at least one of its nearest 
neighbors contains a burning tree. Finally, a burning tree becomes a void in the next timestep. 
Grassberger and Kantz |^ later showed that the model is not critical in the limit p — 0. Later on, 
Drossel and Schwabl ||^ modified the model by introducing a lightning probability / with which a 
green tree catches fire. This version of FFM shows criticality in the sense that anomalous scaling 
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laws are observed provided that {f /p) <^ p ^ ^ / ^ where u' = 0.58 for two-dimensional FFM 
0. In this Section, we compare our model with the FFM introduced by Drossel and Schwabl. 

The parameters of our model should be chosen reasonably so that our model would resemble 
the FFM. First of all, there is no diffusion in the FFM. So we require D -C 1. The propagation of 
fire in the FFM is always allowed and thus we should have Ec2 <^ -Eci- In 2-dimensional FFM, SOC 
behavior is claimed provided that -C p^^ ^ /^^ It can be interpreted as the time 

taken by a forest fire is much shorter than the time taken by the growing of a tree, which is in turn 
much shorter than the time interval of two successive lightning at the same site. In our diffusive 
and dissipative cellular automaton model, the first criterion is satisfied automatically because it 
takes only one timestep to burn down a cluster. For the second criterion, we should set ^ 1 in 
order to reduce the frequency of burning. On the other hand, Ed should not be too large, or else 
the burnable clusters become too large and hence only large events are observed. This is consistent 
with the FFM in which / cannot be too small in comparison to p. After taking all the above criteria 
into consideration, we choose E^i = 6, Ec2 = 1 and D = with a 512 x 512 periodic square lattice 
in our simulation. Of course, other choices of parameters are possible (e.g. E^i = 6, Ec2 = 0.1 and 
D = 10~^) and the results obtained are qualitatively the same. 

In the FFM, a cluster of size s is defined as a group of s neighboring sites with green trees. 
In contrast, the size s of a burnable cluster in our model is defined as the number of sites in that 
burnable cluster. In Fig. ^ we plot the distribution of cluster size N{s) and the root mean quadratic 
radius R{s). We find that 

A^(s) ~s-^ with r = 2.11 ±0.02, (6) 

and 

R{s) ~ s^/^ with fx = 1.87 ± 0.02. (7) 
While the exponent r is similar to the accurate results obtained by Clar et al. H and Grassberger 



10| for the FFM, the exponent /x differs from the numerical value of 1.96 ± 0.01 obtained by Clar 
et al. recently in their FFM simulations 0. 

We then study the model with different Ed and some numerical results are tabulated in Table |. 
The mean density p of sites with energy content greater than the propagation threshold Ec2 can be 
related to the average burning size (S) by 

{S)r^\p.-Pr (8) 

where pc is the mean density at critical point. Using the measured results we have estimated that 

Pc = 0.410 ±0.002, (9) 
7 =2.57±0.01. (10) 

Although the fitted value of 7 is different, pc is the same as that of the FFM within numerical 
accuracy. 

The greatest difference between our model and the FFM is the probability p{s) of burning for 
a cluster of size s. In the FFM, p{s) is proportional to s as all cluster sites are triggering sites in 
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the sense that they are ready to trigger a fire once hghtning occurs. But in our model, only a few 
percent of the sites in a burnable cluster are triggering sites. Fig. |^ shows that the distribution of 
energy dissipation P{S) scales as S""" with an exp{—S/S^) cutoff. We also plot the distribution of 
burning cluster size D{s) and find that within the accuracy of our simulation, P{S) and D{s) share 
the same exponent in the scaling region. Thus we have s ~ 5* and hence 

D{s) oc exp(-s/s^) with a = 0.85 ± 0.01, (11) 

while it is known that a 1.0 for the FFM. Further numerical experiments show that a decreases 
as D increases in our model (see Fig. In conclusion, although our model is similar to that of the 
FFM, they belong to two different universality classes. 

Although a depends on the value of D, we find within numerical errors that the value of a does 
not change with Ed- In Fig. |^ one can see that different E^i affects the range of scaling region but 
not the scaling exponent of P{S). It is easy to realize that for a larger Ed it takes a longer time 
to trigger a burning and allows the growth of larger clusters. In order to prove the fall off of D{s) 
for large s is not due to the finite size of the system, simulations using larger lattice size have been 
done and the same spectra of D{s) (in terms of the number of sites s) are obtained. Furthermore, 
the distribution of time interval between successive burnings P{T) shows an exponential decay as 
discussed in Section II. We conclude in our simulation that no system-wide dissipation occurs. 

In addition to the RMS radius, we define a maximum elongation d{s) which is the maximum 
distance between any two sites in a cluster of size s to investigate the shape of the clusters. We 
found that 

d{s) ~ s^^^' with /i' = 1.82 ± 0.02. (12) 

Fig. tells us that for small s, both R{s) and d{s) deviate from a straight line. This is because 
of the discreteness of the lattice. For larger value of s, however, they follow power laws of similar 
exponents and the ratio d{s)/R{s) increases to a constant value around 2.8 which is very close to 
the value of a 2-dimensional compact circular object (which equals VS). This concludes that small 
burnable clusters are generally irregular and asymmetric objects due to the random fluctuation of 
their formation. As the size of a burnable cluster increases, clusters become more and more regular 
and spherical. The ratio d{s)/R{s) shows that in the average the clusters are more or less compact 
objects for large s. On the other hand, the values of fi and fi' suggest that the clusters may have 
fractal structures with a dimension around 1.9. However, at this stage we cannot determine the 
fractal dimension with accuracy because there is no a general efficient method to estimate it and 
also our simulation is not close enough to the critical point of the system. It is interesting to note 
that when we choose a random site near the center of mass of a large cluster, the average local 
density of the cluster around that site is close to the site percolation threshold 0.59. We beheve 
that there is a close relationship between site percolation and our model (or the FFM) and further 
investigation is needed to draw any conclusion. 

V. CONCLUSIONS AND OUTLOOK 

In this paper, we introduce a cellular automaton model for diffusive and dissipative systems, 
and we have described its behaviors in various different locations in the parameter space. First, we 
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study the effect of diffusion and find tliat P{S) may follow a localized skewed peak, power law, or 
exponential decay. Then we consider the case without diffusion and compare it with the famous 
forest-fire model. 

With a careful choice of parameters and geometry, this model is able to describe qualitatively 
the behavior of a number of physical systems, including the behaviors of type-I X-ray bursters 
(details can be found in Reference 0) and CO2 gas outbursts in some crater lakes (details can be 
found in Reference |^). The following is a brief description of the above two physical systems and 
the reason why our cellular automaton model can be used to describe their behaviors. 

Type-I X-ray bursts are results of thermonuclear explosion on the surface of an accreting neutron 
star. Nuclear fuel, usually making up of hydrogen and helium, is deposited on the neutron star 
surface. Whenever the local density of nuclear fuel is greater than a threshold value, which is 
determined by the specific nuclear reaction involved, a thermonuclear runaway takes place. The 
nuclear reaction can spread to its neighbors if they contain sufficient amount of nuclear fuel. It 
results in a transient flash in the X-ray band, called a type-I X-ray burst. Type-I bursts are different 
from the more frequently repeating type-II bursts, which are believed to be results of accretion disk 
instability [|TT|. It is easy to estimate also that the material diffusion timescale in this problem 



is so long as compared to the typical time interval between two successive bursts. Thus material 
diffusion is not important in this problem. After taken into account for the spatial fluctuation in the 
accretion process, the system can be described using our cellular automaton model on a spherical 
surface with D = 0. In addition, the values of E^i and can be determined in principle once the 
specific nuclear reaction, cooling process, and the typical mass of an accreting blob are given. And 
the burst statistics predicted by our model is consistent with observations 0. 

Another application of the diffusive and dissipative model is the gas outburst statistics in some 
crater lakes. CO2 gas is injected to the bottom of the lake by some natural processes |jl2|. Note 
that the solubility of CO2 in water increases with pressure (and hence the depth of water). So once 
the CO2 concentration in the bottom of the lake becomes supersaturated, gas bubble will form and 
it will drive the water around it to move up causing further degassing. A catastrophic gas outburst 
is resulted. Taking into account the horizontal water flow in the lake, we can describe the outburst 
of a crater lake using our model with a large D, and hence rather regular gas outburst is expected. 

Further investigation of this model in other spatial dimensions, behavior of the model at crit- 
icality as a function of diffusion constant D, and other possible applications of the model will be 
reported in future work 
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APPENDIX: METHOD OF HANDLING DIFFUSION 



The energy difFusion equation in continuous flat two-dimensional space is 

where p is the energy density. In our cellular automaton model, we approximate the above equation 
using the finite difference: 

E(r,t + At) = D^^[E(r',t)-E(r,t)] + E(r,t) for all r (14) 

r' 

where -E(r, t) is the energy content of site r are time t, and the sum is over all the nearest neighbors 
of r. This is possible since the areas of lattice sites are the same. The value of D is coupled with the 
system size and the timestep employed. Eq. ([T^) is a good approximation to the actual diffusion 
process in continuous spacetime if Dtdj IS.I? <^ 1. This criteria is satisfied in all the cases we have 
reported in this paper. In the event that D is too large, an adaptive integration scheme to calculate 
E{r,t + At) is employed. 

The diffusion constant is an invariant under any rescaling of the system. This is true provided 
that the changes in AL and At are correlated. That is. At — > X'^At whenever AL — > XAL. 
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TABLES 



Eel r {E) p 

4 -0.011 ±0.006 53.0 ±0.4 0.3343 ± 0.0005 

5 -0.008 ±0.006 178 ±2 0.363 ± 0.001 

6 +0.012 ±0.006 643 ±6 0.381 ± 0.001 

7 +0.02 ±0.01 2440 ±40 0.393 ± 0.001 



TABLE I. Some numerical results on a 512 x 512 lattice with Ec2 = 1 and D = 0. 



11 



FIGURES 



FIG. 1. A greyscale snapshot of a 512 x 512 lattice with Ed = 6, Ec2 = 1, and D = 0. The higher the 
value of E, the darker the dot is. 

FIG. 2. The variation of the mean dissipation size {S) (solid line) and the correlation coefficient F (dash 
line) against the diffusion constant D. All the simulations are done in a 64 x 64 lattice with Ed = 5.0 and 
Ec2 = 2.0 respectively. The same set of parameters is used in Figs. |3| and ^. 

FIG. 3. Distributions of (a) energy dissipation size S, and (b) time interval between successive energy 
dissipation T for different values of D. 

FIG. 4. Correlation of the time since last energy dissipation T and the size of energy dissipation S for 
different values of D. 

FIG. 5. The distribution of energy dissipation size of a system with D = 5 x 10~^ and Ec2 = 2 (solid 
line). The effect of increasing D to 2 x 10~^ (dotted line) is qualitatively similar to that of reducing Ec2 
to 1 (dashed line). 

FIG. 6. The distribution of (a) cluster size N{s); (b) RMS radius R{s) and maximum elongation d{s). 

FIG. 7. A plot of (a) P{S) against S; and (b) D{s) against s with Ed = 6, Ec2 = 1 and D = in a 
512 X 512 lattice. A scaling region from 10 to 10^ is observed with an exponent of —0.85. 

FIG. 8. Distribution of energy dissipation size P{S) with Ed =7 (solid line), 6(long dash line), 5(dot- 
ted line) and 4 (med dash line) for a 512 x 512 lattice. Ec2 = 1, and D = 0. As since in the S vs. 
Peci{S) / PEcii^O) plot, the scaling region increases with increasing Ed- However, the scaling exponent is 
independent of Ed- 
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